Mound formation and coarsening from a nonlinear instability in surface growth 
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We study a class of one-dimensional, nonequilibrium, conserved growth equations for both non- 
conserved and conserved noise statistics using numerical integration. An atomistic version of these 
growth equations is also studied using stochastic simulation. The models with nonconserved noise 
statistics are found to exhibit mound formation and power-law coarsening with slope selection for a 
range of values of the model parameters. Unlike previously proposed models of mound formation, 
the Ehrlich-Schwoebel step-edge barrier, usually modeled as a linear instability in growth equations, 
is absent in our models. Mound formation in our models occurs due to a nonlinear instability in 
which the height (depth) of spontaneously generated pillars (grooves) increases rapidly if the initial 
height (depth) is sufficiently large. When this instability is controlled by the introduction of an 
infinite number of higher-order gradient nonlinearities, the system exhibits a first-order dynamical 
phase transition from a rough self-affine phase to a mounded one as the value of the parameter that 
measures the effectiveness of control is decreased. We define a new "order parameter" that may be 
used to distinguish between these two phases. In the mounded phase, the system exhibits power-law 
coarsening of the mounds in which a selected slope is retained at all times. The coarsening exponents 
for the continuum equation and the discrete model are found to be different. An explanation of this 
difference is proposed and verified by simulations. In the growth equation with conserved noise, 
we find the curious result that the kinetically rough and mounded phases are both locally stable in 
a region of parameter space. In this region, the initial configuration of the system determines its 
steady-state behavior. 
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I. INTRODUCTION 

The process of growing films by the deposition of atoms 
on a substrate is of considerable experimental and theo- 
retical interest [1]. While there has been a lot of research 
on the process of kinetic roughening [1-3] leading to a 
self-affinc interface profile, there has been much recent 
experimental [4-12] and theoretical [4,6,13-25] interest in 
a different mode of surface growth involving the forma- 
tion of "mounds" which are pyramid-like or "wedding- 
cake-like" structures. The precise experimental condi- 
tions that determine whether the growth morphology 
would be kinetically rough or dominated by mounds are 
presently unclear. However, many experiments show the 
formation of mounds that coarsen (the typical lateral size 
of the mounds increases) with time. During this process, 
the typical slope of the sides of the pyramid-like mounds 
may or may not remain constant. If the slope remains 
constant in time, the system is said to exhibit slope selec- 
tion. As the mounds coarsen, the surface roughness char- 
acterized by the root- mean-square width of the interface 
increases. Eventually, at very long times, the system is 
expected to evolve to a single-mound structure in which 
the mound size is equal to the system size. 

There are obvious differences between the structures 
of kinetically rough and mounded interfaces. In the first 
case, the interface is rough in a self-affine way at length 
scales shorter than a characteristic length ^{t) that ini- 



tially increases with time and eventually saturates at a 
value comparable to the sample size. In the second case, 
the characteristic length is the typical mound size R{t) 
whose time-dependence is qualitatively similar to that 
of (^{t). However, the interface in this case looks well- 
ordered at length scales shorter than R{t). Nevertheless, 
there are certain similarities between the gross features 
of these two kinds of surface growth. First consider the 
simpler situation in which the slope of the sides of the 
mounds remains constant in time. Simple geometry tells 
us that if the system evolves to a single-mound structure 
at long times, then the "roughness exponent" a must 
be equal to unity. Also, the height-difference correla- 
tion function g{r) is expected to be proportional to r for 
r ^ R{t). This is consistent with a = 1. If the mound 
size R{t) increases with time as a power law, R{t) ^ i", 
during coarsening, then the interface width W, which 
is essentially the height of a typical mound, should also 
increase with time as a power law with the same expo- 
nent n. Thus, dynamic scaling with "growth exponent" 
/3 equal to n, and "dynamical exponent" z equal 1/n is 
recovered. If the mound slope s{t) increases with time as 
a power law, s{t) ^ (this is known in the literature as 
steepening) , then one obtains behavior similar to anoma- 
lous dynamical scaling [26] with /3 = n + 6, z ~ 1/n. 

These similarities between the gross scaling properties 
of kinetic roughening with a large value of a and mound 
formation with power-law coarsening make it difficult to 
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experimentally distinguish between these two modes of 
surface growth. This poses a problem in the interpreta- 
tion of experimental results [11,12]. Existing experiments 
on mound formation show a wide variety of behavior. 
Without going into the details of individual experiments, 
we note that some experiments show mound coarsening 
with a time-independent "magic" slope, whereas other 
experiments do not show any slope selection. The de- 
tailed morphology of the mounds varies substantially 
from one experiment to another. The reported values 
of the coarsening exponent n show a large variation in 
the range 0.15-0.4. 

Traditionally, the formation of mounds has been 
attributed to the presence of the so-called Ehrlich- 
Schwoebel (ES) step-edge barrier [27,28] that hinders the 
downward motion of atoms across the edge of a step. This 
step-edge diffusion bias makes it more likely for an atom 
diffusing on a terrace to attach to an ascending step than 
to a descending one. This leads to an effective "uphill" 
surface current [29] that has a destabilizing effect, lead- 
ing to the formation of mounded structures as the atoms 
on upper terraces are prevented by the ES barrier from 
coming down. 

This destabilizing effect is usually represented in con- 
tinuum growth equations by a linear instability. Such 
growth equations usually have a "conserved" form in 
which the time-derivative of the height is assumed to be 
equal to the negative of the divergence of a nonequilib- 
rium surface current j. The effects of an ES barrier are 
modeled in these equations by a term in j that is pro- 
portional to the gradient of the height (for small values 
of the gradient) with a positive proportionality constant. 
Such a term is manifestly unstable, leading to unlimited 
exponential growth of the k 7^ Fourier components of 
the height. This instability has to be controlled by other 
nonlinear terms in the growth equation in order to obtain 
a proper description of the long-time behavior. A number 
of different choices for the nonlinear terms have been re- 
ported in the literature [4,6,13,14,22,23]. If the "ES part" 
of j has one or more stable zeros as a function of the slope 
s, then the slope of the mounds that form as a result of 
the ES instability is expected to stabilize at the corre- 
sponding value(s) of s at long times. The system would 
then exhibit slope selection. If, on the other hand, this 
part of j does not have a stable zero, then the mounds are 
expected to continue to steepen with time. Analytic and 
numerical studies of such continuum growth equations 
have produced a wide variety of results, such as power- 
law coarsening and slope selection with n = 1/4 [13] or 
n ~ 0.17 [6] in two dimensions, power-law coarsening ac- 
companied by a steepening of the mounds [4,14,16], and 
a complex coarsening process [22,23] in which the growth 
of the mound size becomes extremely slow after a char- 
acteristic size is reached. 

There are several atomistic, cellular-automaton-type 
models [19-21] that incorporate the effects of an ES diffu- 
sion barrier. Formation and coarsening of mounds in the 
presence of an ES barrier have also been studied [22,23] 



in a Id model with both discrete and continuum fea- 
tures. We also note that a new mechanism for moiind- 
ing instability has been discovered recently [17,18]. This 
instability, generated by fast diffusion along the edges 
of monatomic steps, leads to the formation of quasi- 
regular shaped mounds in two or higher dimensions. The 
effects of this instability have been studied in simula- 
tions [17,18,20,24,25]. The wide variety of results [17-25] 
obtained from simulations of different models, combined 
with similar variations in the experimental results, have 
made it very difficult to identify the microscopic mecha- 
nism of mound formation in surface growth. 

In this paper, we show that mound formation, slope 
selection and power-law coarsening in a class of one- 
dimensional (Id) continuum growth equations and dis- 
crete atomistic models can occur from a mechanism that 
is radically different from the ones mentioned above. 
Our study is based on the conserved nonlinear growth 
equation proposed by Villain [29] and by Lai and Das 
Sarma [30], and an atomistic version [31] of this equa- 
tion. We have studied the behavior of the continuum 
equation by numerical integration, and the behavior of 
the atomistic model by stochastic simulation. Previous 
work [32-34] on these systems showed that they exhibit 
a nonlinear instability, in which pillars (grooves) with 
height (depth) greater than a critical value continue to 
grow rapidly. This instability can be controlled [32-34] 
by the introduction of an infinite number of higher-order 
gradient nonlinearities. When the parameter that de- 
scribes the effectiveness of control is sufficiently large, 
the controlled models exhibit [32 34] kinetic roughen- 
ing, characterized by usual dynamical scaling with ex- 
ponent values close to those expected from dynamical 
renormalization group calculations [30]. As the value of 
the control parameter is decreased, these models exhibit 
transient multiscaling [32-34] of height fluctuations. For 
yet smaller values of the control parameter, the rapid 
growth of pillars or grooves causes a breakdown of dy- 
namical scaling, with the width versus time plot showing 
a sharp upward deviation [33] from the power-law behav- 
ior found at short times (before the onset of the nonlinear 
instability). 

We report here the results of a detailed study of the be- 
havior of these models in the regime of small values of the 
control parameter where conventional kinetic roughening 
is not observed. We find that in this regime, the inter- 
face self-organizes into a sawtooth-like structiirc with a 
series of triangular, pyramid-like mounds. These mounds 
coarsen in time, with larger mounds growing at the ex- 
pense of smaller ones. In this coarsening regime, a power- 
law dependence of the interface width on time is recov- 
ered. The slope of the mounds remains constant during 
the coarsening process. In section II, the growth equa- 
tion and the atomistic model studied in this work are 
defined and the numerical methods we have used to ana- 
lyze their behavior are described. The basic phenomenol- 
ogy of mound formation and slope selection in these sys- 
tems is described in detail in section III. Specifically, we 
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show that the nonlinear mechanism of mound formation 
in these systems is "generic" in the sense that the qual- 
itative behavior does not depend on the specific form of 
the function used for controlling the instability. In par- 
ticular, we find very similar behavior for two different 
forms of the control function: one used in earlier stud- 
ies [32-34] of these systems, and the other one proposed 
by Politi and Villain [23] from physical considerations. 
Since the linear instability used conventionally to model 
the ES mechanism is explicitly absent in our models, our 
work shows that the presence of step-edge barriers is not 
essential for mound formation. The slope selection found 
in our models is a true example of nonlinear pattern for- 
mation: since the nonequilibrium surface current in our 
models vanishes for all values of constant slope, the se- 
lected value of the slope can not be predicted in any sim- 
ple way. This is in contrast to the behavior of ES-type 
models where slope selection occurs only if the surface 
current vanishes at a specific value of the slope. 

Next, in section IV, wc show that the change in the 
dynamical behavior of the system (from kinetic rough- 
ening to mound formation and coarsening) may be de- 
scribed as a first-order noncqiulibrium phase transition. 
Since the instability in our models is a nonlinear one, 
the fiat interface is locally stable in the absence of noise 
for all values of the model parameters (the strength of 
the nonlinearity and the value of the control parameter). 
The mounded phase corresponds to a different station- 
ary solution of the dynamical equations in the absence of 
noise. We use a linear stability analysis to find the "spin- 
odal boundary" in the two-dimensional parameter space 
across which the mounded stationary solution becomes 
locally unstable. We show that the results of this numer- 
ical stability analysis can also be obtained from simple 
analytic arguments. To obtain the phase boundary in the 
presence of noise, we first define an order parameter that 
is zero in the kinetically rough phase and nonzero in the 
mounded phase. We combine the numerically obtained 
results for this order parameter for different sample sizes 
with finite-size scaling to confirm that this order param- 
eter exhibits the expected behavior in the two phases. 
The phase boimdary that separates the moimdcd phase 
from the kinetically rough one is obtained numerically. 
The phase boundaries for the continuum model with two 
different forms of the control function and the atomistic 
model arc found to be qualitatively similar. 

The results of a detailed study of the process of coars- 
ening of the mounds are reported in section V. Surpris- 
ingly, we find that the coarsening exponents of the con- 
tinuum equation and its atomistic version are different. 
We propose a possible explanation of this result on the 
basis of an analysis of the coarsening process in which the 
problem is mapped to that of of a Brownian walker in an 
attractive force field. In this mapping, the Brownian walk 
is supposed to describe the noise-induced random motion 
of the peak of a mound, and the attractive "force" rep- 
resents the interaction between neighboring mounds that 
leads to coarsening. We show that the numerical results 



obtained for the dynamics of mounds in the atomistic 
model are consistent with this explanation. 

In section VI, we consider the behavior of the con- 
tinuum growth equation for "conserved" noise statistics. 
The nonlinear instability found in the nonconserved case 
is expected to be present in the conserved case also. 
However, there is an important difference between the 
two cases. The nonconserved model exhibits anomalous 
dynamical scaling, so that the typic:al n(;arest-ncighbor 
height difference continues to increase with time, and 
the instability is always reached [33] at sufficiently long 
times, even if the starting configuration is perfectly flat. 
Since the continuum model with conserved noise statis- 
tics exhibits [35] usual dynamic scaling with a < 1, the 
nearest-neighbor height difference is expected to saturate 
at long times if the initial state of the system is flat. Un- 
der these circumstances, the occurrence of the nonlinear 
instability in runs started from flat states would depend 
on the values of the parameters. Specifically, the insta- 
bility may not occur at all if the value of the nonlinear 
coefficient in the growth equation is sufficiently small. At 
the same time, the instability can be initiated by choos- 
ing an initial state with sufficiently high (deep) pillars 
(grooves) . Since mound formation in these models is cru- 
cially dependent on the occurrence of the instability, the 
arguments above suggest that the nature of the long- 
time steady state reached in the conserved model may 
depend on the choice of the initial state. Indeed, we find 
from simulations that in a region of parameter space, the 
mounded and kinetically rough phases are both locally 
stable and the steady state configuration is determined 
by the choice of the initial configuration of the interface. 
These results imply the surprising conclusion that the 
long-time, steady-state morphology of a growing inter- 
face, as well as the dynamics of the process by which the 
steady state is reached may be "history dependent" in 
the sense that the behavior would depend crucially on 
the choice of the initial state. A summary of our find- 
ings and a discussion of the implications of our results 
are provided in Sec. VII. A summary of the basic results 
of our study was reported in a recent Letter [36] . 



II. MODELS AND METHODS 

Conserved growth equations (deterministic part of the 
dynamics having zero time derivative for the k = 
Fourier mode of the height variable) with nonconserved 
noise are generally used [2] to model nonequilibrium sur- 
face growth in molecular beam epitaxy (MBE). The con- 
servation is a consequence of absence of bulk vacancies, 
overhangs and desorption (evaporation of atoms from 
the substrate) under optimum MBE growth conditions. 
Thus, integrating over the whole sample area gives the 
number of particles deposited. This conservation is not 
strictly valid because of "shot noise" fluctuations in the 
beam. The shot noise is modeled by an additive noise 
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term 77(r, t) in the equation of motion of the interface. 
The noise 77 is generally assumed to be delta-correlated 
in both space and time: 



(1) 



where r is a point on a d-dimensional substrate. Thus, a 
conserved growth equation may be written in a form 



(2) 



where h{r,t) is the height at point r at time t, and j is 
the surface current density. The surface current mod- 
els the deterministic dynamics at the growth front. As 
mentioned in section I, the presence of an ES step-edge 
barrier is modeled in continuum equations of the form 
of Eq.(2) by a term in j that is proportional to the slope 
s = V/i, with a positive constant of proportionality. This 
makes the flat surface {h{r) constant for all r) linearly un- 
stable. This instability is controlled by the introduction 
of terms involving higher powers of the local slope s and 
higher-order spatial derivatives of h. 

We consider the conserved growth equation proposed 
by Villain [29] and Lai and Das Sarma [30] for describing 
MBE-type surface growth in the absence of ES step-edge 
barriers. This equation is of the form 

dh'{r,t')/dt' = -v^^h' + A'V2|V/i'|2 + V(r,t'), (3) 

where h'{r,t') represents the height variable at the point 
r at time t'. This equation is believed [2] to provide 
a correct description of the kinetic roughening behavior 
observed in MBE-type experiments [12]. 

In our study, we numerically integrate the Id version 
of Eq.(3) using a simple Euler scheme [33]. Upon choos- 
ing appropriate units of length and time and discretizing 
in space and time, Eq.(3) is written as [33] 

hi{t + At) - hi{t) = AtW^[-'V^hi{t) + A|V/ii(t)|2] 

+ \^tr]i{t), (4) 

where hi{t) represents the dimensionless height variable 
at the lattice point i at dimensionless time t, V and 
are lattice versions of the derivative and Laplacian oper- 
ators, and rii{t) is a random variable with zero average 
and variance equal to unity. Tliesc; equations, with an 
appropriate choice of At, are used to numerically follow 
the time evolution of the interface. In most of our stud- 
ies, we have used the following definitions for the lattice 
derivatives: 



V/ii = {hi+i - hi-i)/2, 
V'^hi = hi+i + hi-i - 2hi. 



(5) 



Wc have checked that the use of more accurate, left-right 
symmetric definitions of the lattice derivatives, involving 
more neighbors to the left and to the right [33], leads to 



results that are very similar to those obtained from cal- 
culations in which these simple definitions arc used. We 
have also checked that the results obtained in the deter- 
ministic limit (?7 = 0) by using a more sophisticated in- 
tegration routine [37] closely match those obtained from 
the Euler method with suflaciently small values of the 
integration time step. 

We have also studied an atomistic version [31] of Eq.(3) 
in which the height variables {hi} arc integers. This 
model is defined by the following deposition rule. First, 
a site (say i) is chosen at random. Then the quantity 



K^{{hJ}) = -V^h, + \\Vh, 



(6) 



is calculated for the site i and all its nearest neighbors. 
Then, a particle is added to the site that has the smallest 
value of K among the site i and its nearest neighbors. In 
the case of a tie for the smallest value, the site i is cho- 
sen if it is involved in the tie; otherwise, one of the sites 
involved in the tie is chosen randomly. The number of 
deposited layers provides a measure of time in this model. 

It was found in earlier studies [32-34] that both these 
models exhibit a nonlinear instability in which isolated 
structures (pillars for A > 0, grooves for A < 0) grow 
rapidly if their height (depth) exceeds a critical value. 
This instability can be controlled [32 -34] by replacing 
[V/ii]^ in Eqns.(4) and (6) by /([V/ii]^) where the non- 
linear function f{x) is defined as 
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(7) 



c > being a control parameter. Wc call the result- 
ing models "model I" and "model II" , respectively. This 
replacement, amounts to the introduction of an infinite 

series of higher-order nonlinear terms. The time evolu- 
tion of the height variables in model I is, thus, given by 

hi{t + At) - hi{t) = AtV^l-V^hiit) 



+ A(1 



-c\Vh,{t)\' 



)/c] + VAtVi{t). (8) 



In model II, the quantity Ki is defined as 

Ki{{hj}) = -V^hi + A(l - e-'=l^^^l')/c. (9) 

While the function f{x) was introduced in the earlier 
work purely for the purpose of controlling the nonlin- 
ear instability, it turns out that the introduction of this 
function in the growth equation is physically meaningful. 
Politi and Villain [23] have shown that the nonequilib- 
rium surface current that leads to the V^jV/i']^ term in 
Eq.(3) should be proportional to V]V/i']^ when |V/i'| is 
small, and should go to zero when JV/i'] is large. The 
introduction of the "control function" /([V/ij]^) satisfies 
this physical requirement. We have also carried out stud- 
ies of a slightly different model (which we call "model 
lA") in which the function f{x) is assumed to have a 
form suggested by Politi and Villain: 
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where c is, as before, a positive control parameter. This 
function has the same asymptotic behavior as that of 
the function defined in Eq.(7). As we shall show later, 
the results obtained from calculations in which these two 
different forms of /(.x) arc used arc qualitatively very 
similar. In fact, we expect that the qualitative behavior 
of these models would be the same for any monotonic 
function f{x) that satisfies the following requirements: 
(i) f{x) must be proportional to x in the small- a; limit; 
and (ii) it must saturate to a constant value as a; oo. 

Wc have carried out extensive simulations of both 
these models for different system sizes. The results 
reported here have been obtained for systems of sizes 
40 < L < 1000. There is no significant dependence of 
the results on L. The time step used in most of our stud- 
ies of models I and lA is At ~ 0.01. We have checked 
that very similar results are obtained for smaller values of 
At. We used both unbounded (Gaussian) and bounded 
distributions for the random variables rji in our simula- 
tions of models I and lA, with no significant difference in 
the results. For computational convenience, a bounded 
distribution (uniform between +\/3 and — \/3) was used 
in most of our calculations. Unless otherwise stated, the 
results described in the following sections were obtained 
using periodic boundary conditions. The efi^ects of using 
other boundary conditions will be discussed in the next 
section. 



III. MOUND FORMATION AND SLOPE 
SELECTION 

It has been demonstrated earlier [32,33] that if the 
control parameter c is sufficiently large, then the nonlin- 
ear instability is completely suppressed and the models 
exhibit the usual dynamical scaling behavior with the 
expected [30] exponent values, /3 ~ 1/3, 2: ~ 3, and 
a = /32 ~ 1. This behavior for model I is illustrated 
by the solid line in Fig.l, which shows a plot of the 
width 14^ as a function of time t for parameter values 
A=4.0 and c = 0.06. As the value of c is decreased with 
A held constant, the instability makes its appearance: 
the height ho of an isolated pillar (for A > 0) increases 
in time if hmin{^,c) < ho < hjnax{^,c). The value of 
hmin is nearly independent of c, while hmax increases 
as c is decreased [33]. If c is sufficiently large, h„iax is 
small and the instability does not affect the scaling be- 
havior of global quantities such as W, although transient 
multiscaling at length scales shorter than the correlation 
length ^ ~ t^/^ may be found [32,33] if c is not very 
large. As c is decreased further, hmax becomes large, 
and when isolated pillars with ho > hmin ai"© created at 
an initially flat interface through fluctuations, the rapid 
growth of such pillars to height hmax leads to a sharp 
upward departure from the power-law scaling of W with 



time t. The time at which this departure occurs varies 
from run to run. This behavior for model I with A = 4.0 
and c — 0.02 is shown by the dash-dotted line in Fig.l. 

This instability leads to the formation of a large num- 
ber of randomly distributed pillars of height close to 
hmax- As the system evolves in time, the interface self- 
organizes to form triangular mounds of a fixed slope near 
these pillars. These mounds then coarsen in time, with 
large mounds growing larger at the expense of small ones. 
In this coarsening regime, a power-law growth of W in 
time is recovered. The slope of the sides of the triangu- 
lar mounds remains constant during this process. Finally, 
the system reaches a steady state with one peak and one 
trough (if periodic boundary conditions are used) and 
remains in this state for longer times. The interface pro- 
files in the kinetically rough phase (obtained for relatively 
large values of c) and the mounded phase (obtained for 
small c) are qualitatively different. This difference is il- 
lustrated in Fig. 2 that shows typical interface profiles in 
the two different phases. This figure also shows a typical 
interface profile for model I A in the mounded regime, il- 
lustrating the fact that the precise choice of the control 
function f{x) is not crucial for the formation of mounds. 
The evolution of the interface structure in the mounded 
regime of model I is illustrated in Fig. 3 which shows 
the interface profiles obtained in a typical L = 200 run 
starting from a flat initial state at three different times: 
t = 200 (before the onset of the instability), t = 4000 (af- 
ter the onset of the instability, in the coarsening regime), 
and t = 128000 (in the final steady state). This figure 
also shows the steady-state interface profile of a L = 500 
sample with the same parameters, to illustrate that the 
results do not depend on the sample size. 

Very similar behavior is found for model II. Since the 
heights in this atomistic model can increase only by dis- 
crete amounts in each unit of discrete time, the increase 
of W at the onset of the instability is less rapid here than 
in the continuum models I and lA. Nevertheless, the oc- 
currence of the instability for small values of c shows 
up as a sharp upward deviation of the W versus t plot 
from the initial power-law behavior with f3 ~ 1/3. This 
is illustrated by the dash-dotted line in Fig. 4, obtained 
from simulations of model II with A = 2.0, c = 0.005. 
This behavior is to be contrasted with that for A = 2.0, 
c = 0.015, shown by the full line in Fig. 4, where the 
nonlinear instability is absent. The difference between 
the surface morphologies in the two regimes of model II 
is illustrated in Fig. 5. The kinetically rough, self-affine 
morphology obtained for c = 0.02 is clearly different from 
the mounded profile found for c = 0.005. The time evolu- 
tion of the interface in the mounded regime of this model 
is illustrated in Fig. 6. The general behavior is clearly 
similar to that found for models I and lA. This figure 
also shows a properly scaled plot of the interface profile 
of a i = 500 sample with the same parameters at a time 
in the coarsening regime. It is clear from this plot that 
the nature of the interface and the value of the selected 
slope do not depend on the sample size. 
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The occiirrcnce of a peak and a symmetrically placed 
trough in the steady-state profiles shown in Figs 3 and 6 
is a consequence of using periodic boundary conditions. 
The deterministic part of the growth equation of Eq.(8) 
strictly conserves the average height if periodic boundary 
conditions are used. So, the average height remains very 
close to zero if the initial state is flat, as in most of our 
simulations. The steady-state profile must have at least 
one peak and one trough in order to satisfy this require- 
ment. Also, it is easy to show that if the slopes of the 
"uphill" and "downhill" parts of the steady-state profile 
are the same in magnitude (this is true for our mod- 
els), then the two extrema must be separated by ~ L/2. 
Therefore, it is clear that the steady state obtained in 
simulations with periodic boundary conditions must have 
a peak and a symmetrically placed trough separated by 
distance ~ L/2. 

To check whether the basic phenomenology described 
above depends on the choice of the boundary condition, 
we have carried out test simulations using two other 
boundary conditions: "fixed" boundary condition, in 
which the height variable to the left of i = 1, and to the 
right oi i = L are pinned to zero at all times; and " zero 
flux" boundary condition with vanishing flrst and third 
derivatives of the height at the two ends of the sample. 
For these boundary conditions, the deterministic part of 
the growth equation does not strictly conserve the av- 
erage height. As a result, the symmetry between the 
mound and the trough, found in the long-time steady 
state obtained for periodic boundary condition, is not 
present if one of the other boundary conditions is used. 
In particular, it is possible to stabilize a single mound or 
a single trough in the steady state for the other boundary 
conditions. Since the heights at the two ends must be the 
same for flxed boundary condition, the two extrema in a 
configuration with one mound and one trough must be 
separated by ~ i/2, as shown in Fig. 7. The two extrema 
would not be separated by ~ L/2 for the zero-flux bound- 
ary condition. These effects of boundary conditions are 
illustrated in Fig. 7 which shows profiles in the mounded 
regime obtained for the three different boundary condi- 
tions mentioned above. It is clear from the results shown 
in this figure that the basic phenomenology, i.e. the for- 
mation and coarsening of mounds and slope selection, is 
not affected by the choice of boundary conditions. In par- 
ticular, the values of the selected slope and the heights 
of the pillars at the top of a mound and the bottom of 
a trough remain unchanged when boundary conditions 
other than periodic are used. 

The selection of a "magic slope" during the coars- 
ening process is clearly seen in the plots of Fig. 3 and 
Fig. 6. More quantitatively, the probability distribution 
of the magnitude of the nearest-neighbor height differ- 
ences Si = \hi+i — hi\ is found to exhibit a pronounced 
peak at the selected value of the slope, and the position 
of this peak docs not change during the coarsening pro- 
cess. Fig. 8 shows a comparison of the distribution of the 
magnitude of the nearest-neighbor height difference for 



model I in the mounded and kinctically rough phases. 
A bimodal distribution is seen for the mounded phase, 
the two peaks corresponding to the values of the selected 
slope and the height of the pillars at the top and bottom 
of the pyramids. The kinetically rough phase, on the 
other hand, exhibits a distribution peaked at zero. Fig. 9 
shows the values of the selected slope at different times 
in the coarsening regime of model I. The constancy of the 
slope is clearly seen in this plot. All these features remain 
true for the discrete model. Plots of the distribution P{s) 
at two different times in the coarsening regime of model 
II are shown in Fig. 10. The peak position shows a small 
shift in the positive direction as t is increased, but this 
shift is small compared to the width of the distributions, 
indicating near constancy of the selected slope. The value 
of the selected slope depends on the parameters A and c. 
This is discussed in the next section. 



IV. DYNAMICAL PHASE TRANSITION 

The instability that leads to mound formation in our 

models is a nonlinear one, so that the perfectly flat state 
of the interface is a locally stable steady-state solution of 
the zero-noise growth equation for all parameter values. 
When the instability is absent (e.g. for large values of 
the control parameter c), this "fixed-point" solution of 
the noise-free equation is transformed to the kinetically 
rough steady state in the presence of noise. The mounded 
steady state obtained for small values of c must corre- 
spond to a different fixed point of the zero-noise growth 
equation. Such nontrivial fixed-point solutions may be 
obtained from the following simple calculation. 

The profile near the top {i = ig) of a triangular 
mound may be approximated a,s hi^ = xq + Xi, hi^^+j = 
2^0 — (IjI — l)a^2, where x\ is the height of the pillar at 
the top of the mound and X2 is the selected slope. This 
profile would not change under the dynamics of Eq.(8) 
with no noise if the following conditions are satisfied: 

V2/i,„-A(l-e-^l^''*ol')/c 
= V%„±i-A(l-e-^l^''*o±ir)/c 
= V2/iio±2 - A(l - e-^l^'*^o±2|=)/c. (11) 

These conditions lead to the following pair of non-linear 
equations for the variables xi and X2 used to parametrize 
the profile near the top of a mound: 

2xi-A[l-e-'=^']/c=0, 
3a;i -X2- A[l - e-^^^^^+^^^'/^J/c = 0. (12) 

These equations admit a non-trivial solution for suffi- 
ciently small c, and the resulting values of xi and x^ are 
found to be quite close to the results obtained from nu- 
merical integration. A similar analysis for the profile near 
the bottom of a trough (this amounts to replacing X2 by 
—X2 in Eq.(12)) yields slightly different values for x\ and 
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X2. The full stable profile (a fixed point of the dynamics 
without noise) with one peak and one trough may be ob- 
tained numerically by calculating the values of \hi\ for 
which (/j, the term multiplying At in the right-hand side 
of Eq.(8), is zero for all %. The fixed-point values of {/i,} 
satisfy the following equations: 

5i = \]\-\]-^hi + A(l - e-'=l'^'''l')/c] = for all i. (13) 

A numerical solution of these coupled nonlinear equations 
shows that the small mismatch between the values of xi 
near the top and the bottom is accommodated by creat- 
ing a few ripples near the top. The numerically obtained 
fixed-point profile for a i = 500 system with A = 4.0, 
c = 0.02 is shown in Fig.ll, along with a typical steady- 
state profile for the same system. The two profiles are 
found to be nearly identical, indicating that the mounded 
steady state in the presence of noise corresponds to this 
fixed-point solution of the noiseless discretized growth 
equation. 

Fixed-point solutions of the continuum equation, 
Eq.(3), with i/ = 1 and |Vft,p replaced by /(|Vft,|2) where 
/(x) has the form shown in Eq.(lO) may also be obtained 
by a semi-analytical approach following Racz tt al. [38]. 
We consider stationary solutions of the continuum equa- 
tion that satisfy the following first-order differential equa- 
tion with appropriate boundary conditions: 



where s{x) = dh{x)/dx is the local slope of the interface 
and A is a constant that must be positive in order to 
obtain a solution that resembles a triangular mound. At 
large distances from the peak of the mound, the slope 
s would be constant, so that ds{x)/dx would vanish, 
whereas the second term would give a positive contri- 
bution if A is positive. At the peak of the profile, the sec- 
ond term would be zero because s is zero, but ds{x)/dx 
would be negative, making the left-hand side of Eq.(14) 
positive. While a closed-form solution of this differen- 
tial equation cannot be obtained, the value of s{x) at 
any point x may be calculated with any desired degree of 
accuracy by numerically solving a simple algebraic equa- 
tion. The height profile is then obtained by integrating 
s{x) with appropriate boundary conditions. In our calcu- 
lation, we used the procedure of Racz et. al. [38] to take 
into account periodic boundary conditions. In Fig. 11, 
we have shown a typical steady-state profile of a L = 200 
sample of model lA with A = 4.0 and c = 0.01, and 
a fixed-point solution of the corresponding continuum 
equation. The value of the constant A in Eq.(14) was 
chosen to yield the same slope as that of the steady-state 
profile of the discrete model. These results show that the 
steady-state properties for the two forms of f{x) are very 
similar, and the continuum equation admits stationary 
solutions that are very similar to those of the discretized 
models. 



The local stability of the mounded fixed point may be 
determined from a calculation of the eigenvalues of the 
stability matrix, Afy = dgi/dhj, evaluated at the fixed 
point. We find that the largest eigenvalue of this matrix 
(disregarding the trivial zero eigenvalue associated with 
an uniform displacement of the interface, hi ^ hi + S for 
all i) crosses zero at c = ci(A) (see Fig. (12)), signaling 
an instability of the mounded profile. The structure of 
Eq.(8) implies that ci(A) cx A^. Thus, for < c < ci(A), 
the dynamics of Eq.(8) without noise admits two locally 
stable invariant profiles: a trivial, flat profile with hi the 
same for all i, and a non-trivial one with one mound and 
one trough. Depending on the initial state, the noise- 
less dynamics takes the system to one of these two fixed 
points. For example, an initial state with one pillar on 
a fiat background is driven by the noiseless dynamics to 
the flat fixed point if the height of the pillar is smaller 
than a critical value, and to the mounded one otherwise. 

The "relevant" perturbation that makes the mounded 
fixed point unstable at c = ci is a uniform vertical rela- 
tive displacement of the segment of the interface between 
the peak and the trough of the fixed-point profile. This 
can be seen by numerically evaluating the right eigen- 
vector corresponding to the eigenvalue of the stability 
matrix that crosses zero at c = ci . This is demonstrated 
in the inset of Fig. 12. Also, examination of the time evo- 
lution of the mounded structure for values of c slightly 
higher than ci shows that the instability of the struc- 
ture first appears at the bottom of the trough. Taking 
cue from these observations, the value ci can be obtained 
from a simple calculation. We consider the profile near 
the bottom of a trough at « = iq. As discussed above, 
the profile near zq may be parametrized as hi^ = xq+xi, 
^ia+j = '^0 + (b'l — 1)2^2, and the values of xi and X2 
may be obtained by solving a pair of nonlinear equa- 
tions, Eq.(12) with X2 replaced by — .T2. Wo now con- 
sider a perturbation of this profile, in which the heights 
on one side of io are all increased by a small amount 5 
(i.e. hig+j = a;o + (j - ^)x2 + S, hi^^j = xq + (j - l)x2 
with j > 0), and use Eq.(8) to calculate how 5 changes 
with time, assuming its value to be small. The require- 
ment that 5 must decrease with time for the fixed-point 
structure to be locally stable leads to the following equa- 
tion for the value of c at which the structure becomes 
unstable: 

^{xi - a;2)e-'=(^^-^=)'/4 = 1, (15) 

By substituting the numerically obtained values of Xi and 
X2 in this equation, the critical value, ci (A) , of the param- 
eter c is obtained as a function of A. The values obtained 
this way are in good agreement with those obtained from 
our full numerical calculation of the eigenvalues of the 
stability matrix. The "spinodal" lines (i.e. the lines in 
the c — A plane beyond which the mounded fixed point 
is unstable) for models I and lA are shown in Fig. 13. 
Both these lines have the expected form, ci(A) oc A^. It 
would be interesting to carry out a similar stability anal- 
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ysis for the mounded stationary profile (see Fig. 11) of the 
continuum equation corresponding to model lA. Such a 
calculation would have to be performed without discretiz- 
ing space if we want to address the question of whether 
the behavior of the truly continuum equation is similar 
to that of the discretized versions considered here. We 
have not succeeded in carrying out such a calculation: 
since the mounded stationary profiles for the continuum 
equation are obtained from a numerical calculation, it 
would be extremely difficult, if not impossible, to carry 
out a linear stability analysis for such stationary solutions 
without discretizing space. 

In the presence of the noise, the perfectly flat fixed 
point transforms to the kinetically rough steady state, 
and the non-trivial fixed point evolves to the mounded 
steady state shown in Fig.ll. A dynamical phase tran- 
sition at c = 02(A) < ci(A) separates these two kinds 
of steady states. To calculate C2(A), we start a system 
at the mounded fixed point and follow its evolution ac- 
cording to Eq.(8) for a long time (typically t = 10**) to 
check whether it reaches a kinetically rough steady state. 
By repeating this procedure many times, the probability, 
Pi (A, c), of a transition to a kinetically rough state is ob- 
tained. For fixed A, Pi increases rapidly from to 1 as 
c is increased above a critical value. Typical results for 
Pi as a function of c for model I with A = 4.0 are shown 
in the inset of Fig. 13. The value of c at which Pi = 0.5 
provides an estimate of C2 . Another estimate is obtained 
from a similar calculation of P2(A, c), the probability that 
a flat initial state evolves to a mounded steady state. As 
expected, P2 increases sharply from to 1 as c is de- 
creased (see inset of Fig. 13), and the value of c at which 
this probability is 0.5 is slightly lower than the value at 
which Pi = 0.5. This difference reflects flnite-timc hys- 
teresis effects. The value of C2 is taken to be the average 
of these two estimates, and the difli'erence between the 
two estimates provides a measure of the uncertainty in 
the determination of C2. The phase boundary obtained 
this way is shown in Fig. 13, along with the results for 
C2(A) obtained for the discrete model II from a similar 
analysis. 

The general behavior found for all the models as the 
parameters A and c are varied is qualitatively very simi- 
lar to that in equilibrium first order phase transitions of 
two- and three-dimensional systems as the temperature 
and other parameters, such as the magnetic field in spin 
systems, are varied. To take a standard example of an 
equilibrium first order transition, we consider a system 
with a scalar order-parameter field ip{r), described by a 
Ginzburg-Landau free energy functional [39] that has a 
cubic term: 



/(m) = -arn^ — -bm^ + -um^. 

^ O 41 



(17) 



I 



dr 



(16) 



where b and u are positive constants, o = oo(T— Tq) with 

flg > 0, and T is the temperature. Considering uniform 
states, '0(r) = m, the free energy per unit volume may 
be written as 



It is easy to show that for To < T < = Tq + b'^ / (Aaou) , 
the function f{m) has two local minima, one at m = 0, 
and the other at a positive value of m. These two minima 
represent the two phases of the system. This system ex- 
hibits a first order equilibrium phase transition from the 
disordered phase (m = 0) to an ordered phase with pos- 
itive m as the temperature is decreased. The transition 
temperature lies between Tq and Tg. The temperature 
Ts at which the minimum corresponding to the ordered 
phase disappears is called the "spinodal" temperature for 
the ordered phase. The spinodal temperature for the dis- 
ordered phase is Tq. 

Now consider the dynamics of this system according 
to the following time-dependent Ginzburg-Landau equa- 
tion [39]: 



dip{r, t) 
dt 



(18) 



where F is a kinetic coefficient and r] represents Gaus- 
sian delta-correlated noise whose variance is related to 
F and the temperature T via the fiuctuation-dissipation 
theorem [39] . In the absence of noise, this equation con- 
verges to local minima of the functional F. So, the 
noiseless dynamics exhibits two locally stable fixed points 
for To < T < Ts, corresponding to the two minima of 
f{m) that represent the disordered and uniformly or- 
dered states. This is analogous to the two locally stable 
fixed points of our nonequilibrium dynamical systems for 
c < Ci(A). If we identify the flat and mounded fixed 
points as the "disordered" and "ordered" states, respec- 
tively, and the control parameter c to play the role of the 
temperature T, then the noiseless dynamics of our mod- 
els would look similar to that of Eq.(18) for ij = 0, with 
ci playing the role of the spinodal temperature T^ of the 
equilibrium problem. 

In the presence of noise, the system described by 
Eq.(18) exhibits a first-order phase transition at Tc that 
lies between T, and Tq: the system selects one of the 
phases corresponding to the two fixed points of the noise- 
less dynamics, except at T^. where both phase coexist. 
The local stability of the mean-field ordered and disor- 
dered states in a small temperature-interval around Tc 
is manifested in the dynamics as finite-time hysteresis 
effects. The behavior we find for out nonequilibrium 
dynamical models is qualitatively similar: the system 
selects the steady state corresponding to the mounded 
("ordered") fixed point of the noiseless dynamics as the 
control parameter c (analogous to the temperature T of 
the equilibrium system) is decreased below C2 which is 
smaller than the "spinodal" value ci. The growth mod- 
els do not exhibit a "spinodal" point for the kinetically 
rough ("disordered") phase: the fiat fixed point of the 
noiseless dynamics is locally stable for all positive values 
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of the control parameter c. If this analogy with equi- 
librium first order transition is correct, then our mod- 
els should show hysteresis and coexistence of kinetically 
rough and mounded morphologies for values of c near 
C2(A). As mentioned above, we do find hysteresis (see 
inset of Fig. 13) in finite-time simulations with values of c 
near C2. Evidence for two-phase coexistence is presented 
in Fig. 14, where a snapshot of the interface profile for a 
L = 500 sample of model I with A = 4.0, c = 0.42 is 
shown. This value of c is very close to the critical value 
C2 for A = 4.0 (see inset of Fig. 13). This plot clearly illus- 
trates the simultaneous presence of mounded and rough 
morphologies in the interface profile. 

The results described above suggest that our growth 
models exhibit a first-order dynamical phase transition 
at c = C2(A). To make this conclusion more concrete, 
we need to define an order parameter, analogous to the 
quantity m in the equilibrium problem discussed above, 
that is zero in the kinetically rough phase, and jumps 
to a non-zero value as the system undergoes a transition 
to the mounded phase at c = C2. The identification of 
such an order parameter would also be useful for dis- 
tinguishing between these two different kinds of growth 
in experiments - as mentioned in the Introduction, it is 
difficult to experimentally differentiate between kinetic 
roughening and mound formation with coarsening from 
measurements of the usual bulk properties of the inter- 
face. A clear distinction between the two morphologies 
may be obtained from measurements of the average num- 
ber of extrema of the height profile [40] . The steady-state 
profile in the mound-formation regime exhibits two ex- 
trema for all values of the system size L. In contrast, 
the number of extrema in the steady state in the kinetic 
roughening regime increases with i as a power law [40] 
we find that for values of c for which the system is kinet- 
ically rough, e.g. for A = 4.0, c = 0.05 for model I, the 
average number of extrema in the steady state is propor- 
tional to L'' with 5 ~ 0.83. This observation allows us 
to define an "order parameter" that is zero in the large- 
c, kinetic roughening regime and finite in the small-c, 
mound-formation regime. Let (Ji be an Ising-like variable, 
equal to the sign of the slope of the interface at site i. 
An extremum in the height profile then corresponds to a 
"domain wall" in the configuration of the {ui} variables. 
Since there are two domain walls separated by ^ L/2 
in the steady state in the mound-formation regime, the 
quantity 

m=^|(X:a,-e2-^/^)|, (19) 

where (. . .) represents a time-average in the steady state, 
would be finite in the L oo limit. On the other hand, 
m would go to zero for large L in the kinetically rough 
regime because the number of domains in the steady- 
state profile would increase with L. We find numeri- 
cally that in the kinetically rough phase, m ~ with 



7 ~ 0.2. The finite-size scaling data for the order param- 
eter m for models I and II for both faceted and kineti- 
cally rough phases is shown in Fig. 15. It is seen that mL 
varies linearly with the system size L in the mounded 
phase, whereas mL ~ L^~^ with 7 ~ 0.2 for model I and 
7 ~ 0.15 for model II in the the kinetically rough phase. 
So, in the L ^ 00 limit, the order parameter would jump 
from zero to a value close to unity as c is decreased be- 
low C2(A). This is exactly the behavior expected at a 
first-order phase transition. 

The occurrence of a first-order phase transition in our 
Id models with short-range interactions may appear sur- 
prising - it is well-known [39] that Id systems with short- 
range interactions do not exhibit any equilibrium thermo- 
dynamic transition at a non-zero temperature. The situ- 
ation is, however, different for nonequilibrium phase tran- 
sitions: In contrast to equilibrium systems, a first-order 
phase transition may occur in one- dimensional nonequi- 
librium systems with short-range interactions. Several 
such transitions have been well documented in the lit- 
erature [41]. So, there is no reason to a priori rule out 
the occurrence of a true first-order transition in our Id 
nonequilibrium systems. As discussed above, our numer- 
ical results strongly suggest the existence of a true phase 
transition. However, since all our results are based on 
finite-time simulations of finite-size systems, we can not 
claim to have established rigorously the occurrence of a 
true phase transition in our models. The crucial question 
in this context is whether the order parameter m would 
be nonzero in the mounded phase in the L oo limit if 
the time-average in Eri.(19) is performed over arbitrarily 
long times. Since the steady-state profile in this phase 
has a single mound and a single trough (this is clear from 
our simulations) , the only way in which m can go to zero 
is through strong "phase fluctuations" corresponding to 
lateral shifts of the positions of the peak and the trough. 
We do not find any evidence for such strong phase fluc- 
tuations. We have calculated the time autocorrelation 
function of the phase of the order parameter for small 
samples over times of the order of 10^ and found that it 
remains nearly constant at a value close to unity over the 
entire range of time. So, if such phase fluctuation even- 
tually make the order parameter zero for all values of c, 
then this must happen over astronomically long times. 
Our flnite-time simulations can not, of course, rule out 
this possibility. 



V. COARSENING OF MOUNDS 

During the late-stage evolution of the interface, the 
mounds coarsen with time, increasing the typical size 
of the triangular pyramidal structures. The process of 
coarsening occurs by larger mounds growing larger at the 
expense of the smaller ones while always retaining their 
"magic" slope. Snapshots of the system in the coarsen- 
ing regime are shown in Figs 16 and 17 for model I and 
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model II, respectively. The constancy of the slope during 
the coarsening process is clearly seen in these figures. As 
discussed in the Introduction, the constancy of the slope 
implies that if the typical lateral size of a mound increases 
in time as a power law with exponent n {R{t) oc t"), then 
the width of the interface would also increase in time as 
a power law with the same exponent {W{t) oc with 
(3 = n). Therefore, the value of the coarsening expo- 
nent n may be obtained by measuring the width as a 
function of time in the coarsening regime. In Fig. 18, we 
show a plot of the width as a function of time for model 
I with A = 4.0, c = 0.02. It is clear from the plot that 
the time-dependence of the width is well-described by a 
power law with (3 = n = 0.34 ± 0.01. A similar plot for 
the discrete model II with A = 2.0, c = 0.005, shown 
in Fig. 19, also shows a power-law growth of the width 
in the long-time regime, but the value of the coarsen- 
ing exponent obtained from a power-law fit to the data 
is (3 = n = 0.50 ±0.01, which is clearly different for 
the value obtained for model I. This is a surprising re- 
sult: model II was originally defined [31] with the spe- 
cific purpose of obtaining an atomistic realization of the 
continuum growth equation of Eq.(3), and earlier stud- 
ies [31-33] have shown that the dynamical scaling be- 
havior of this model in the kinetic roughening regime is 
the same as that of model I. Also, we have found in the 
present study that the dynamical phase transition in this 
model has the same character as that in model I. So, the 
difference in the values of the coarsening exponents for 
these two models is unexpected. As noted earlier, there 
is some evidence suggesting that the typical slope of the 
mounds in model II increases very slowly with time (see 
Fig. 10). However, this "steepening" , if it actually occurs, 
is too slow to account for the large difference between the 
values of the coarsening exponents for models I and II. 

In order to understand these numerical results, we first 
address the question of why the mounds coarsen with 
time. This problem has certain similarities with domain 
growth in spin systems [42]. Using the Ising variables 
{cTj} defined in the preceding section, each height profile 
can be mapped to a configuration of Ising spins. The 
coarsening of mounds then corresponds to a growth of 
the typical size of domains of these Ising spins. There 
is, however, an important difference between the coars- 
ening of mounds in our models and the usual domain 
growth problem [42] for Ising spins. Domain growth in 
spin systems is the process through which the system ap- 
proaches equilibrium from an out-of-equilibrium initial 
state. The dynamics of this process may be understood 
in terms of arguments based on considerations of the free 
energy (at finite temperatures) or energy (at zero tem- 
perature). Such arguments do not apply to our nonequi- 
librium growth models. The reason for the coarsening of 
mounds in our models must be sought in the relative sta- 
bility of different structures under the assumed dynamics 
and the effects of noise. 

As discussed in the preceding section, the fixed point 
of Eq.(8) with one mound and one trough is locally sta- 



ble for c < ci(A). Since structures with several mounds 
and troughs approach this steady-state structure through 
the coarsening process, it is reasonable to expect that 
fixed points of the noiseless equation with more than one 
mounds and troughs would be unstable. We have numer- 
ically obtained fixed points of Eq.(8) with two mounds 
and troughs for different values of the sample size and the 
separation between the peaks of the two mounds. The 
slope of the mounds at these fixed points is found to be 
the same as that in the fixed point with one mound and 
one trough. We find that the stability matrix for such 
fixed points always has a real, positive eigenvalue, indi- 
cating that the structure is unstable and would evolve to 
the stable configuration with one mound and one trough. 
The magnitude of the positive eigenvalue of the stability 
matrix for two-mounded fixed points depends on the sam- 
ple size, the separation between the peaks of the mounds 
and the relative heights of the mounds in a complicated 
way. We have not been able to extract any systematic 
quantitative information from these dependences. We 
find a qualitative trend indicating that the magnitude 
of the positive eigenvalue decreases as the separation be- 
tween the peaks of the two mounds is increased. Since the 
time scale of the development of the instability of two- 
mounded structures is given by the inverse of the positive 
eigenvalue, this result is consistent with the expectation 
that the time required for two mounds to coalesce should 
increase with the separation between the mounds. 

These results suggest that the coarsening of the 
mounds in model I reflects the instability of structures 
with multiple mounds and troughs. If this is true, then 
coarsening of mounds should be observed in this model 
even when the noise term in Eq.(8) is absent. To check 
this, we have carried out numerical studies of coarsening 
in the noiseless version of Eq.(8). In these studies, the 
time evolution of an initial configuration with a pillar of 
height ho > hmini^, c) at the central site of an otherwise 
flat interface is followed numerically in the presence of 
noise until the instability caused by the presence of the 
pillar is well developed. The profiles obtained for differ- 
ent realizations of the noise used in the initial time evolu- 
tion arc then used as initial configurations for coarsening 
runs without noise. The dotted line in Fig. 18 shows the 
width versus time data obtained from this calculation. 
The coarsening exponent in the absence of noise is found 
to be the same {n ~ 1/3) as that of the noisy system, 
indicating that the coarsening in this model is driven by 
processes associated with the deterministic part of the 
growth equation. 

We have examined the details of the process by which 
two mounds coalesce to form a single one. The different 
steps in this process are illustrated in the snapshots of in- 
terface profiles shown in Fig. 17 where one can sec how the 
two mounds near the center come together to form a sin- 
gle one as time progresses. First, the separation between 
the peaks of the mounds decreases with time. When this 
distance becomes sufficiently small, the "V" -shaped seg- 
ment that separates the peaks of the mounds "melts" to 
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form a rough region with many spikes. This region of 
the interface then self-organizes to become the top part 
of a mound. Although the data shown in this figure were 
obtained for model II, it also represents quite closely the 
process of coalescence of mounds in models I and I A. An 
estimate of the time scale associated with the second part 
of this process, during which the "melted" region of the 
interface transforms into the top part of a mound, may 
be obtained in the following way. The absolute value 
of the closest-to-zero eigenvalue of the stability matrix 
for the single-mounded fixed point provides an estimate 
of the time scale over which configurations close to the 
fixed point evolve to the fixed point itself. In the inset of 
Fig. 18, we have shown the dependence of the magnitude 
Km of the closest-to-zero eigenvalue for A = 4.0, c = 0.02 
on the system size L. The eigenvalue scales with the sys- 
tem size as L~^, indicating that the time scale for the 
decay of fluctuations with length scale L is proportional 
to L^. This is consistent with the observed value of the 
coarsening exponent, n ~ 1/3, which indicates that the 
time scale t{x) for the coalescence of mounds separated 
by distance x is proportional to a;^/" ~ x^. 

We have found very similar behavior for the closcst- 
to-zero eigenvalue of the stability matrix for the single- 
mounded fixed point of model lA in which a different 
form of the control function is used. Coarsening data for 
this model are shown in Fig. 18. In this model, there is 
a long time interval between the onset of the instability 
and the beginning of power-law coarsening. During this 
time interval, the interface segments near the tall pillars 
formed at the instability organize themselves into trian- 
gular mounds. This process produces a plateau in the 
width versus time plot. Eventually, however, power-law 
coarsening with n ~ 1/3 is recovered, as shown by the 
dashed line in Fig. 18. Since the onset of power-law coars- 
ening in this model occurs at very late times, we could 
not get coarsening data for this model over a very wide 
time interval. Consequently, the calculated value of the 
coarsening exponent for this model is less accurate. How- 
ever, our results show quite clearly that the "universal" 
features of the coarsening dynamics of models I and lA 
are the same. 

Going back to the discrete model II, we first examined 
its coarsening dynamics in the absence of noise. The 
noiseless limit of this model is not well-defined in the 
sense that there is no explicit noise term that can be 
turned off. The stochasticity in this model arises from 
two sources: first, the randomness associated with the se- 
lection of the deposition site i (the quantity Ki{{hj}) de- 
fined in Eq.(9) is calculated at this site and at its nearest- 
neighbor sites); and second, the randomness in the selec- 
tion of one of the two neighbors of site i in case of a tie in 
their values of Ki. In order to make the dynamics deter- 
ministic, we employ a parallel update scheme in which all 
the lattice sites, i = 1, - ■ ■ L, are updated simultaneously 
instead of sequentially. This eliminates the stochasticity 
arising from the choice of the sequence in which the sites 
are updated. The randomness associated with the selec- 



tion of a neighbor in case of a tie is eliminated by choosing 
the right neighbor if the serial number of the occurrence 
of a tie, measured from the beginning of the simulation, 
is even, and the left neighbor if the serial number is odd. 
With these modifications of the update rules, the system 
evolves in a perfectly deterministic way. To study coars- 
ening in this deterministic version of model II, we prepare 
an initial structure with two identical mounds separated 
by distance xq. The slope of these mounds is chosen to 
be equal to the "selected" value found in simulations of 
the original model. We then study the time evolution 
of this structure according to the parallel dynamics de- 
fined above, monitoring how the distance x between the 
peaks of the mounds changes with time t. We find that 
the the value of x increases initially, in order to accom- 
modate a slightly higher value of the selected slope for 
the parallel dynamics. After reaching a maximum value 
that is ~ 20% higher than xq, the distance x decreases 
with time, indicating that the noiseless dynamics leads 
to coarsening. Eventually, the two mounds coalesce into 
a single one and the system remains in the state with 
one mound and one trough at later times. The slope of 
the mounds remains constant during the coarsening pro- 
cess. Assuming power-law coarsening with exponent n, 
the separation x at time t is expected to have the form 

x{t)^{xl/"-Ctr, (20) 

where C is a constant. This form implies that the time 
t{xo) required for the coalescence of two mounds sepa- 
rated by distance xq is proportional to xj^". In Fig. 20, 
we have shown the time dependence of x{t) for three dif- 
ferent initial separations, xq = 80,90 and 100, and fits 
of the data to the form of Eq.(20), yielding the result 
1/n = 2.9 ±0.1. Only the data for x at times larger 
than the time at which it returns to xq after the initial 
increase are shown in the figure. From these observa- 
tions, we conclude that the coarsening exponent in the 
zero-noise limit of model II is the same (n = 1/3) as that 
found for the two versions of the continuum model. The 
observation that the coarsening exponent for the noisy 
version of model II is different from 1/3 then indicates 
that the effect of noise in the discrete model is qualita- 
tively different from that in the continuum models. We 
discuss below a possible explanation of this behavior. 

The fact that noiseless versions of all three models ex- 
hibit the same value of the coarsening exponent (n = 1/3) 
suggests that the coarsening is driven by an effective 
attractive interaction between the peaks of neighboring 
mounds. The observed value of n suggests [43] that the 
this attractive interaction is proportional to 1/x^ where 
X is the separation between the mound tips. This inter- 
action would lead to the observed result, t{x) oc x'^, in 
the noiseless limit if the rate of change of x with t is as- 
sumed to be proportional to the attractive force ("over- 
damped limit"). The presence of noise in the original 
growth model leads to a noise term in the equation of 
motion of the variable x, but the nature of this noise term 
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is not clear. Since the observed coarsening dynamics in 
the noisy model II {n ~ 0.5) suggests a similarity with 
random walks, we propose that the effective dynamics of 
the variable x is governed by the kinetic equation 

^ = -C/x'+v{t), (21) 

where 77 is a Gaussian, delta-correlated noise with zero 
mean and variance equal to 2D. In this phenomenolog- 
ical description, the coarsening of a two-mounded struc- 
ture in model II is mapped to a Brownian walk of a par- 
ticle in an attractive potential field with an absorbing 
wall at the origin, such that the particle cannot escape 
once it arrives at the origin. The absorption of a par- 
ticle at the origin corresponds to the coalescence of two 
mounds in the original height picture. Thus, in this re- 
duced model, the quantity of interest is the typical first 
passage time t (i.e. the time taken by a particle to reach 
the origin) as a function of xq, the initial distance of the 
particle from the origin. In the noiseless limit {D = 0), 
T is equal to Xq/(3C), and in the purely Brownian walk 
limit (C — 0), the typical value of t should be of the 
order of Xq/D. Therefore, for sufficiently large values 
of xq, random- walk behavior characterized by n = 1/2 
is expected. However, for relatively small values of xq, 
the behavior should be dominated by the attractive in- 
teraction. Therefore, we expect that the dynamics de- 
scribed by Eq.(21) with nonzero C and D would exhibit 
a crossover from a noise dominated regime to an inter- 
action dominated regime as the value of xq is decreased. 
This crossover is expected to occur near xq ~ Xc ^ C/D, 
for which the values of r obtained from the two individ- 
ual terms in Eq.(21) become comparable. We, therefore, 
propose a scaling form for the dependence of r on a;o: 

where the scaling function F[z) has the following asymp- 
totic dependence on its argument z: F{z) = 1 for 2; <C 1 
and F{z) oc I/2; for 2 » 1. Our numerical study of the 
reduced model of Eq.(21) confirms the validity of this 
scaling ansatz. Note that for any nonzero value of D, 
r(xo) would be proportional to Xq, implying n = 1/2, for 
sufficiently large values of xq. 

To test the validity of this reduced description of 
the coarsening dynamics of the original model II with 
stochastic sequential updates, we have simulated the evo- 
lution of a two-mounded structure in this model. The 
two-mounded structure used in these simulations is iden- 
tical to that used in the study of coarsening in the model 
with parallel updates. In these simulations also, the aver- 
age value of x exhibits a small initial increase followed by 
a steady decrease. The initial growth of {x'^{t)) — {x{t))^ 
with time is found to be linear, indicating the presence 
of a random additive noise in the effective cqiiation of 
motion of x. Fig. 21 shows the time dependence of {x{t)) 
obtained from simulations of L = 500 samples of model 



II with A = 2.0, c = 0.005, and a;o = 100. In this plot, 
the origin of time has been shifted to the point where 
(x) returns to the initial value xq — 100 after the small 
initial increase, and 10^ units of simulation "time" (num- 
ber of deposited layers) is taken to be the unit of t. The 
number of independent runs used in the calculation of 
averages is 800. The observed dependence of {x{t)) on 
t can be described reasonably well by the reduced equa- 
tion of Eq.(21) for appropriate choice of the values of 
the parameters C and D. As shown in Fig. 21, the {x{t)) 
calculated numerically from Eq.(21) with C = 285.0 and 
2D = 0.3 provides a good fit to the data obtained from 
simulations of the growth model. Due to the limited 
range of the simulation data, the values of C and D can 
not be determined very accurately from such fits: values 
of C in the range 250 — 300 and values of D in the range 
0.1 — 0.5 (larger values of C have to be combined with 
smaller values of D) provide reasonable descriptions of 
the simulation data. For such values of C and D, and 
xo ~ L where L = 10^ is the sample size used in the 
calculation of the coarsening exponent, Dxq/C is of or- 
der unity, indicating that the effects of the noise term in 
Eq.(21) should be observed in the simulation data. We, 
therefore, conclude that the presence of an additive ran- 
dom noise term in the effective equation of motion for 
the separation between the peaks of neighboring mounds 
in model II is a plausible explanation for the observed 
value of the coarsening exponent, n = 1/2. 

In view of this conclusion, it is interesting to enquire 
why the coarsening exponent n for models I and lA has 
the value 1/3 characteristic of dynamics governed by the 
deterministic interaction between mound tips. We can 
not provide a conclusive answer to this question. One 
possibility is that the additive random noise in the origi- 
nal growth equations for these models does not translate 
into a similar noise in the effective equation of motion 
for the separation of mound tips. A second possibility 
is that the equation of motion for the separation x for 
these models also has the form of Eq.(21), but the rel- 
ative strength of the noise is much smaller, so that the 
crossover value Xc is much larger than the typical sample 
sizes used in our simulations. Under these circumstances, 
the dynamics of x would be governed by the interaction 
and r would be proportional to Xq, giving the value 1/3 
for the coarsening exponent n. If the second explana- 
tion is correct, then one should observe a crossover from 
n = l/3ton = 1/2 in models I and I A as the sample 
size L is increased. We do not find any evidence for such 
a crossover in our simulations. 

In passing, we note that the dynamics of the slope 
variables {s,} is strictly conserved in the models studied 
here if periodic or fixed boundary conditions are used. 
Also, the deterministic part of the growth equations con- 
serves the integrated height. The "magnetization" of the 
Ising variables {ai} representing the signs of {sj} is not 
strictly conserved: it is conserved only in a statistical 
sense. However, unlike other well-known examples [42] 
of systems with conserved dynamics, we obtain in model 
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II a coarsening exponent that is different from the ex- 
pected Lifshitz-Slyozov value, n ~ 1/3. Since the height 
variables in model II are integers, there are some sites 
for which hi = hi+i. The assignment of the value of the 
Ising variable at at such sites is clearly ambiguous. It 
may be more appropriate to use a three-state variable, 
taking the values ±1 and 0, to describe the coarsening 
behavior of this model. This problem does not arise in 
models I and lA for which the height variables are real 
numbers because the likelihood of two neighboring height 
variables to be exactly the same is vanishingly small. 

VI. MODEL I WITH CONSERVED NOISE 

As discussed in section IV, the properties of the 
mounded phase of model I are determined to a large ex- 
tent by the mounded fixed point of the deterministic part 
of the equations of motion of the height variables. The 
presence of noise changes the critical value of the con- 
trol parameter c from ci to C2 < Ci, but does not affect 
strongly the properties of the mounded steady state of 
the system (sec, for example. Fig. 11). Therefore, we ex- 
pect that the properties of the mounded phase would not 
change drastically if the statistics of the noise is altered. 
On the other hand, it is well-known [1,2] that the expo- 
nents that describe the scaling behavior in the kinetically 
rough phase depend strongly on the nature of the noise. 
In particular, the exponents for conserved noise are ex- 
pected to be quite different from those describing the 
behavior for nonconserved noise. The occurrence of the 
nonlinear instability that leads to the mounded phase in 
our models is contingent upon the spontaneous forma- 
tion of pillars of height ho > hmin{^,c), if the initial 
state of the system is completely flat. The probability 
of formation of such pillars depends crucially [33] on the 
values of the roughening exponents which, in turn, are 
strongly dependent on the nature of the noise. There- 
fore, we expect that the nature of the noise may be very 
important in determining whether the instability leading 
to mound formation actually occurs in samples with flat 
initial states. 

We have investigated this issue in detail by carrying 
out simulations of a version of Model I in which the noise 
is conserved [35] . The equations of motion for the height 
variables in this model have the form of Eq.(8), with the 
noise terms {r]i{t)} having the properties 

Mt)) = 0, MtHit')) = -^^kA,t', (23) 

where 6ij = 1 if i = j and zero otherwise. This model 
is expected to exhibit kinetic roughening with exponents 
(3 ~ 1/11, a — 1/3, and z ~ 11/3 in one dimension [35]. 
Since the value of a for this model is less than unity, it 
exhibits conventional dynamical scaling with the typical 
value of the nearest-neighbor height difference saturat- 
ing at a constant value at long times. The value of this 
constant is expected to increase [33] as the strength A of 



the nonlinearity is increased. As discussed in detail in 
Rcf. [33], the nonlinear instability that leads to mound 
formation is expected to occur in the time evolution of 
such models from a perfectly flat initial state only if the 
value of A is sufficiently large to allow the spontaneous, 
noise-induced formation of pillars of height greater than 
hmini^, c). So, if the value of A is sufficiently small, then 
the model with conserved noise, evolving from a flat ini- 
tial state, would not exhibit the mounding transition. On 
the other hand, if the instability is induced in the model 
by starting the time evolution from a state in which there 
is at least one pillar with height greater than h„iin, then it 
is expected to evolve to the mounded state if the value of 
c is sufficiently small to make the mounded state stable. 
So, the long-time steady state of the conserved model is 
expected to exhibit an interesting dependence on the ini- 
tial state: if A is sufficiently small (so that pillars with 
height greater than hmin are not spontaneously generated 
in the time evolution of the interface from a flat initial 
state), and the value of c sufficiently small (so that the 
mounded state is stable in the presence of noise), then 
the steady state would be kinetically rough if the initial 
state is sufficiently smooth, and mounded if the initial 
state contains pillars of height greater than hmin- This 
"bistability" does not occur for the nonconserved model 
I because the nearest-neighbor height difference in this 
model continues to increase with time, so that the insta- 
bility always occurs at sufficiently long times [33]. 

Our simulations of the model with conserved noise 

show the bistable behavior discussed above in a large 
region of the c — A plane. We find that in this model, the 
height of a pillar on an otherwise flat interface increases 
in time if its initial value ho is larger than hmin — 20/ A 
(the dependence of hmin on c is weak) . This dependence 
of hmin on A is very similar to that [33] found for model 
I with nonconserved noise. We also find that the typ- 
ical values of the nearest-neighbor height difference do 
not continue to increase with time in this model. Con- 
sequently, if A is sufficiently small, pillars with height 
greater than hmin are not generated, and the system ex- 
hibits conventional kinetic roughening with exponent val- 
ues close to the expected ones [35]. On the other hand, 
if the time evolution of the same system is started from 
a state with a pillar of height greater than hmim then it 
evolves to a mounded state very similar to the one found 
in the model with nonconserved noise if the value of c 
is sufficiently small. The two steady states obtained for 
the conserved model with A = 4.0, c = 0.01 are shown 
in Fig. 22. The long-time state obtained in a run start- 
ing from a flat configuration is kinetically rough, whereas 
the state obtained in a run in which the nonlinear insta- 
bility is initially seeded in the form of a single pillar of 
height ho = 1000 at the central site is mounded with a 
well-defined slope, as in model I with nonconserved noise. 
The difference between the two profiles, obtained for the 
same parameter values for two different initial states, is 
quite striking. 
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Since the steady state in the conserved model depends 
on the initial condition, it is not possible to draw a con- 
ventional phase diagram for this model in the c — A plane: 
the transition lines are different for different initial con- 
ditions. In Fig. 23, we have shown three transition lines 
for this model in the c — A plane. The line drawn through 
the circles (line 1) is obtained from simulations in which 
the system is started from a flat initial state. If A is 
small, then the steady state reached in such runs is ki- 
netically rough for all values of c. As A is increased 
above a "critical value" Ac — 5.3, pillars with height 
greater than hmin sltc spontaneously generated during 
the time evolution of the system and it exhibits a tran- 
sition to the mounded state if the value of c is not very 
large. The circles represent the values of c for which 
50% of the runs show transitions to the mounded state. 
The line through the diamonds (line 2) corresponds to 
50% probability of transition to the mounded state from 
an initial state with a pillar of height ho = 1000 on an 
otherwise flat interface. The probability of reaching a 
mounded steady state in such runs decreases from unity 
as the value of c is increased, and falls below 50% as 
line 2 is crossed from below. For large A, lines 1 and 2 
merge together. This is expected: the probability of oc- 
currence of a mounded steady state should not depend 
on how the pillars that initiate the nonlinear instability 
are generated. The third line (the one passing through 
the squares) represents 50% probability of transition to 
the kinetically rough state from a mounded initial state 
(the fixed point of the noiseless equations of motion with 
one mound and one trough). This line reflects the noise- 
induced instability of the mounded steady state for rel- 
atively large values of c. The differences between lines 2 
and 3 are due to flnite-time hysteresis effects similar to 
those discussed in section IV in the context of determin- 
ing the critical value C2(A) of the control parameter c for 
model I with nonconserved noise. 

The interesting region in the "phase diagram" of Fig. 23 
is the area enclosed by lines 1 and 2 and the c = line. 
For parameter values in this region, the system exhibits 
bistable behavior, as discussed above. This bistability is 
unexpected in the sense that in most studies of nonequi- 
librium surface growth, it is implicitly assumed that the 
long-time steady state of the system docs not depend 
on the choice of the initial state. So, it is important 
to examine whether the dependence of the steady state 
on the initial condition in the conserved model reflects 
a very long (but flnite) transition time from one of the 
two apparent steady states to the other one. We have 
addressed this question by carrying out long {t of the or- 
der of lO'') simulations of small samples with parameter 
values in the middle of the "bistable region" for flat and 
mounded initial states. We did not find any evidence for 
transitions from one steady state to the other one in such 
simulations. Of course, we can not rule out the possibil- 
ity that such transitions would occur over much longer 
time scales. 



VII. SUMMARY AND DISCUSSIONS 



To summarize, we have shown from numerical simula- 
tions that a class of Id surface growth models exhibits 
mound formation and power-law coarsening of mounds 
with slope selection as a result of a nonlinear instabil- 
ity that is controlled by the introduction of an infinite 
series of terms with higher-order gradient nonlinearities. 
The models considered here arc discrctized versions of a 
well-known continuum growth equation and an atomistic 
model originally formulated for providing a discrete real- 
ization of the growth equation. We have shown that these 
models exhibit a dynamical phase transition between a 
kinetically rough phase and a mounded phase as a pa- 
rameter that measures the effectiveness of controlling the 
instability is varied. We have defined an order parameter 
for this first-order transition and used finite-size scaling 
to demonstrate how the sample-size dependence of this 
order parameter provides a clear distinction between the 
rough and mounded phases. We have also mapped out 
the phase boundary that separates the two phases in a 
two-dimensional parameter space. 

We would like to emphasize that the ES mechanism, 
commonly believed to be responsible for mound forma- 
tion in surface growth, is not present in our models. Our 
models exhibit a nonlinear instability, instead of the lin- 
ear instability used conventionally to represent the ef- 
fects of ES barriers. The mechanism of mound forma- 
tion in our models is also different from a recently dis- 
covered [17,18] one involving fast edge diffusion, which 
occurs in two or higher dimensions. The slope selection 
found in our models is a rare example of pattern for- 
mation from a nonlinear instability. This is clearly dif- 
ferent from slope selection in ES-type models in which 
the mounds maintain a constant slope during coarsening 
only if the nonequilibrium surface current vanishes at a 
particular value of the slope. The selected slope in such 
models is simply the slope at which the current is zero. 
The behavior of our models is more complex: in these 
models, the siuface current is zero for all values of con- 
stant slope, and the selected value of the slope is obtained 
from a nonlinear mechanism of pattern selection. 

Our studies bring out two other unexpected results. 
We find that the coarsening behavior of an atomistic 
model (model II) specifically designed to provide a dis- 
crete realization of the growth equation that leads to 
model I is different from that of model I: the exponents 
that describe the power-law coarsening are different in 
the two models. We show that this difference may arise 
from a difference in the nature of the effective noise that 
enters the equation of motion for the separation between 
neighboring mounds in the two cases. The second sur- 
prising result is that the numerically obtained long-time 
behavior of model I with conserved noise in a region of 
parameter space depends crucially on the initial condi- 
tions: the system reaches a mounded or kinetically rough 
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steady state depending on whether or not the initial state 
is sufficiently rough. To our knowledge, this is the first 
example of "nonergodic" behavior in nonequilibrium sur- 
face growth. 

The behavior found in our Id models may be relevant 
to experimental studies of the roughening of steps on a 
vicinal surface. As noted earlier, the form of the con- 
trol function /(x) used in model I A is physically rea- 
sonable. However, since very little is known about the 
values of the model parameters appropriate for exper- 
imentally studied systems, we are unable to determine 
whether the mechanism of mound formation found in 
our study would be operative under experimentally re- 
alizable conditions. The nonlinear instability found in 
our Id models is also present [34] in the experimentally 
interesting two-dimensional version of the growth equa- 
tion of Eq.(3). However, it is not clear whether this in- 
stability, when controlled in a manner similar to that in 
our Id models, would lead to the formation of mounds 
in two dimensions. This question is currently under in- 
vestigation. Since the growth equation of Eq.(3) exhibits 
conventional dynamic scaling in the kinetic roughening 
regime in two dimensions, the nonlinear instability would 
not occur in runs with flat initial states if the value of 
A is small. Therefore, the behavior in two dimensions is 
expected to be similar to that of our Id model I with 
conserved noise: the nature of the long-time steady state 
may depend crucially on initial conditions in a region of 
parameter space. Such nonergodic behavior, if found in 
two dimensions, would have interesting implications for 
the growth of films on patterned substrates. 

All the results described in this paper have been ob- 
tained from numerical studies of models that are dis- 
crete in both space and time. It is interesting to enquire 
whether the truly continuum growth equation of Eq.(3) 
exhibits similar behavior. This question acquires special 
significance in view of studies [33,44] that have shown 
that discretization may drastically change the behavior 
of nonlinear growth equations similar to Eq.(3). Since 
the interesting behavior found in our discretized models 
arises from the nonlinear instability found earlier [32,33], 
the question that we have to address is whether a sim- 
ilar instability is present in the truly continuum growth 
equation. This question was addressed in some detail in 
Ref. [33] where it was shown that the nonlinear instabil- 
ity is not an artifact of discretization of time or the use 
of the simple Euler scheme for integrating the the growth 
equation. In the present study, we have found additional 
evidence (see section HI) that supports this claim. We 
should also point out that the atomistic model II, for 
which the question of inaccuracies arising from time in- 
tegration does not arise, exhibits very similar behavior, 
suggesting that the behavior found in models I and lA 
is not an artifact of the time discretization used in the 
numerical integration. 

The occurrence of the nonlinear instability does de- 
pend on the way space is discretized (i.e. how the lattice 
derivatives are defined). In earlier work [32,33] as well 



as in the present study, the lattice derivatives are de- 
fined in a left-right symmetric way. We have found that 
the instability actually becomes stronger if the number 
of neighbors used in the definition of the lattice deriva- 
tives is increased. This result suggests that the instability 
is also present in the continuum equation. It has been 
found by Putkaradze et al. [45] that the instability does 
not occur if the lattice derivatives are defined in a dif- 
ferent way in which cither left- or right-discretization of 
the nonlinear term is used, depending on the sign of the 
local slope of the interface. There is no reason to believe 
that this definition is "better" or "more physical" than 
the symmetric definitions used in our work. The only rig- 
orous result we arc aware of for the behavior of Eq.(3) is 
derived in Ref. [45] where it is shown that the solutions of 
the noiseless equation are bounded for sufficiently smooth 
initial conditions. This result, however, does not answer 
the question of whether the instability occurs in the con- 
tinuum equation. As discussed in Ref. [33] , the nonlinear 
instability of Eq.(4), signalled by a rapid initial growth 
of the height (depth) of isolated pillars (grooves), may 
not lead to a true divergence of the height variables. The 
results reported in the present work would remain valid 
as long as high pillars or deep grooves are formed by the 
instability - the occurrence of a true divergence is not 
necessary. 

In the present work, we have shown (see section IV) 
that the continuum equation with f{x) defined in Eq.(lO) 
does admit stationary solutions that exhibit all the rel- 
evant features of stationary solutions of the discretized 
equation. This result provides additional support to the 
contention that the behaviors of the continuum and dis- 
cretized systems are qualitatively similar. We should, 
however, mention that these stationary solutions of the 
continuum equation do not pick out a selected slope of 
the interface: profiles similar to those shown in Fig. 11 
may be obtained for different values of the parameter 
A in Eq.(14) that determines the slope of the triangu- 
lar mound. Slope selection in the continuum equation 
may occur as a consequence of the requirement of local 
stability of such stationary solutions. As mentioned in 
section IV, wc have not attempted a linear stability anal- 
ysis of such numerically obtained stationary solutions of 
the continuum equation because doing such a calculation 
without discrctizing space would be extremely difficult. 
Further investigation of this question would be useful. 

Finally, we would like to emphasize that the discrete 
models studied here would continue to be valid models for 
describing nonequilibrium surface growth even if the be- 
havior of the truly continuum growth equation of Eq. (3) 
turns out to be different from that found here. These 
models may be looked upon as ones in which continu- 
ous (in models I and lA) or discrete (in model II) height 
variables defined on a discrete lattice evolve in continu- 
ous or discrete time. These models have all the correct 
symmetries and conservation laws of the physical prob- 
lem, and they exhibit, for different values of the control 
parameter c, both the phenomena of kinetic roughening 
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and mound formation found in experiments. There is no 
compelling reason to consider the continuum equation 
to be more "correct" or "physical" than these models. 
Epitaxial growth is intrinsically a discrete process at the 
molecular level and a continuum description is an ap- 
proximation that may not be valid in some situations. 

Prom a different perspective, the nonequilibrium first- 
order phase transition found in our models is interest- 
ing, especially because it occurs in Id systems with short 
range interactions. Such phase transitions have been 
found earlier in several Id "particle hopping" models [41]. 
It would be interesting to explore possible connections 
between such models and the Id growth models studied 
here. 
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FIG. 1. Double-log plots of the interface width W as a func- 
tion of time t for model I with A = 4.0, c — 0.02 (dash-dotted 
line), and for A = 4.0, c = 0.06 (full line). These data are for 
L = 500 samples, averaged over 40 independent runs starting 
from flat initial states. Plots have been shifted vertically for 
clarity. 



FIG. 3. The interface profile at three different times {t — 
200, 4000, and 128000) in a run starting from a flat state for 
a L = 200 sample of model I with A=4.0 and c=0.02. The 
dashed line shows the profile for a L = 500 sample with the 
same parameters at t = 1.28 x lO'^, with both axes scaled by 
2.5. 
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FIG. 2. Configuration snapshots at t=W^, for model I with 

A=4.0, c=0.06 (top panel), model I with A=4.0, c=0.02 (mid- 
dle panel), and model lA with A=4.0, and c=0.01 (bottom 
panel). 



FIG. 4. Double-log plots of the interface width as a func- 
tion of time t for model II with A=2.0, c=0.005 (dash-dotted 
line), and for A=2.0, c=0.015 (full line). These data are for 
L — 500 samples, averaged over 40 independent runs starting 
from flat states. Plots have been shifted vertically for clarity. 



17 




FIG. 
A=2.0, 
panel). 



5. Configuration snapshots at t—10 
c=0.015 (top panel), and A=2.0, 



for model II with 
c=0.005 (bottom 
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FIG. 6. The interface profile at three times (t = 10^, 10®, 
and 3 x 10®) in a run starting from a flat state for a L = 200 
sample of model II with A=2.0 and c=0.005. The profile of 
a L = 500 sample with the same parameters at t = 3 x 10®, 
with both axes scaled by 2.5, is shown by the dashed line. 



FIG. 7. Height profiles for model I (A = 4.0, c = 0.02) at 
time t — 1.28 x 10^ for periodic boundary conditions (full 
line), fixed boundary conditions (dashed line), and zero-flux 
boundary conditions (dash-dotted line). 
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FIG. 8. Distribution of the magnitude of the near- 
est-neighbor height difference s for model I with A=4.0 and 
c=0.02 (left panel), showing the bimodal nature of the dis- 
tribution, characteristic of a mounded phase with slope selec- 
tion. The distribution for A=4.0, c=0.05 (right panel) does 
not show this behavior. 
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FIG. 9. Average slope of the mounds as a function of time 
for model I in the mounded phase (A — 4.0, c — 0.02) during 
the coarsening process, i=8000 to t=1.28 x 10^. The data are 
for L=500 samples averaged over 40 runs. 
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FIG. 11. Fixed point profile for a, L — 500 sample of model 
I with A=4.0 and c=0.02 (full line), compared with a steady 
state profile (dash-dotted line) for the same parameter values. 
The dashed line shows a steady state profile of a L = 200 sam- 
ple of model lA with A = 4.0, c = 0.01 (both axes scaled by 
2.5), and the dotted line shows an invariant solution of the 
corresponding continuum equation. 
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FIG. 10. Distribution of the magnitude of the near- 
est-neighbor height difference s for model II with A=2.0 and 
c=0.005, at two different times, t=10^ (full line) and t=10^ 
(dashed line). 
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FIG. 12. The dependence of k, the closest-to-zero eigen- 
value of the stability matrix for the mounded fixed point 
of model I with A=4.0, L=500, on the control parameter c. 
The inset shows the right eigenvector p, corresponding to this 
eigenvalue near the point where k crosses zero. 
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FIG. 13. Critical values of the control parameter c as func- 
tions of A: ci of model I (circles), C2 of model I (triangles), 
C2 of model II (diamonds), and ci of model lA (squares). In- 
set: The probabilities Pi (circles) and P2 (triangles) defined 
in text, as functions of c for model I with A = 4.0, L = 200. 
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FIG. 14. Two phase coexistence near the phase transition 
in model I. The plot shows an interface profile of a L = 500 
sample with A=4.0, c=0.42. 



FIG. 15. Finite-size scaling for the order parameter m. 
The left panel shows double-log plots of mL as a function 
of the sample size L for model I in the mounded phase 
(A=4.0, c=0.02, circles) with slope 1.00 ± 0.01 and in the ki- 
netically rough phase (A=4.0, c=0.05, diamonds) with slope 
0.81 ± 0.02. The right panel shows similar plots for model 
II in the mounded phase (A=2.0, c=0.005, circles) with 
slope 1.00 ± 0.01 and in the kinetically rough phase (A=2.0, 
c=0.015, diamonds) 0.88 ± 0.02. The straight lines are the 
best power-law fits to the data. 
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FIG. 16. Snapshots (t = 2 x 10*, 6 x 10"*, 10^ 1.4 x 10^ 
and 1.28 x lO'^) of the profile of a L = 500 sample of model 

I (profiles at different times have been shifted in the vertical 
direction for clarity) with A=4.0, c=0.02 in the coarsening 
regime. 
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FIG. 17. Snapshots {t = 1100000, 1200000, 1250000, 
1301000, 1450000, 1606000, 1660000, 1670000, 1680000, 
1700000) of the profile of a L = 500 sample of model II 
(profiles at different times have been shifted in the vertical 
direction for clarity) with A=2.0, c=0.005 in the coarsening 
regime. 




L ^ In L 7 

10 20 

FIG. 18. Double-log plots of the interface width W as a 

function of time t for model I with A=4.0 and c=0.02, with 
noise (dash-dotted line), and without noise (dotted line) av- 
eraged over 40 runs for L — 1000 samples. The dashed line 
shows W versus t data for model lA with A=4.0, c=0.01, 
averaged over 60 runs for L = 500 samples. The solid line 
represents power- law behavior with exponent n = 1/3. Inset: 
Finite-size scaling data for the inverse of the closest-to-zero 
eigenvalue of the stability matrix for the mounded fixed point 
of model I with A = 4.0, c = 0.02. 



FIG. 19. Double-log plot of the interface width W as a 
function of time t for model II with A=2.0 and c=0.005 
(solid line) for L=1000 samples averaged over 40 runs. The 
dash-dotted and dashed lines represent power-law behavior 
with exponent n = 1/3 and n = 1/2, respectively. 



X 




FIG. 20. Peak separation function of time t for a 

two-mounded structure for model II with parallel updates (see 
text). The data shown are for A=2.0, c=0.005, L=500. The 
initial value of the separation is xo=80 (squares), a;o=90 (di- 
amonds), and a;o = 100 (circles). The solid lines represent the 
fits described in the text. 
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FIG. 21. The solid line shows the average value of the sepa- 
ration x{t) between mound tips (see text) as a function of time 
t for a two-mounded structure for model II. The data showu 
are for A=2.0, c=0.005, L=5QQ, a;o=100, averaged over 800 
runs. The dashed line shows {x{t)) calculated for the reduced 
model of Eq.(21) with C = 285.0 and D = 0.15. 
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FIG. 22. Profiles at time t=10^ for model I with conserved 

noise (A=4.0, c=0.02, L — 500), for a flat initial configuration 
(top panel) , and an initial configuration with a pillar of height 
ho = 1000 (bottom panel). 




FIG. 23. Phase diagram for model I with conserved noise. 
The plots show 50% stability lines (sec text) for a fiat initial 
state (circles), an initial state with a pillar of height ho = 1000 
(diamonds) and an initial state identical to the mounded fixed 
point of the noiseless equations of motion (squares). The data 
were obtained from 100 t = 10* runs for L = 200 samples. 
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